Load data

Intensity and rain data is handled here.

Load Rain Data

Rain is cut to the narrow area surrounding the wetlands as plotted in the time series.

## although coordinates are longitude/latitude, st_intersects assumes that they are planar

As stars, Assign date dimension

Set polygonID. Data is clipped accordingly

Load reference stars objects

all <- loadPolygon(study_area[1,])
# don't know why but can only give vector here
# all_bigger <- loadPolygon(study_area[2,])
# lose the extra dimension
# all_bigger <- all_bigger[,,,1]

Create DF with VV measures & classes

end <- length(water_shape$type)
for (k in 1:end) {
  if(water_shape[k,]$type != "mixed") {
    poly <- loadPolygon(water_shape[k,])
    # set VV, timestep
    poly <- poly[1,,,1]
    poly.df <- as.data.frame(poly)
    poly.df$type <- water_shape[k,]$type
    poly.df <- poly.df[complete.cases(poly.df),]
    if(k == 1) {
      val <- poly.df[,4:5]
    }
    else {
      val <- rbind(val, poly.df[,4:5])
    }
  }
}
# boxplot
ggplot(val, aes(x=type, y=VV)) +
  geom_boxplot() +
  geom_hline(yintercept = -15.34187)

Implement SVM-like Algorithm and calculate Threshold

# thresh <- findThreshold(val)
thresh <- list("threshold" = -14.77885)
thresh
## $threshold
## [1] -14.77885

Plot AOI with calculated Threshold

# plot overall region
image(all[1,,,1], col = c("white", "black", "green"), breaks = c(-30, thresh$threshold, -8.833493, 20))

Export as Raster

Time Series: Rain, Scaled Plot and Threshold

Rain is summed up over the plotted area. Absolute values are given in mm and relative values are given in respect to the yearly sum.

plotTS(all, thresh$threshold, -8.833493, 1, 31)
## Warning: Removed 11 rows containing missing values (position_stack).

Plots

colo = viridisLite::inferno(31)
breaks = seq(0, 31, 1)
image(class.sum, col = colo, breaks = breaks)
image(class.sum[study_area[1,]], col = colo, breaks = breaks)
plot(class.sum[moni[12,]])
plot(class.sum[shape[124,]])

SVM

Prepare new data without street class..

# prepare data
val.svm$type <- as.factor(val.svm$type)
attach(val.svm)
x <- subset(val.svm, select=-type)
y <- type
# make model
svmmod <- svm(x,y)
summary(svmmod)
## 
## Call:
## svm.default(x = x, y = y)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1 
## 
## Number of Support Vectors:  24
## 
##  ( 12 12 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  field water
# make prediction, confusion matrix
pred <- predict(svmmod,x)
table(pred, y)
##        y
## pred    field water
##   field   246     2
##   water     0   826
# threshold is
svmmod$x.scale$`scaled:center`
## [1] -15.34187
# plot classified with threshold calculated by SVM
plot(all[1,,,1], col = c("white", "black"), breaks = c(-35, svmmod$x.scale$`scaled:center`, 10))
# vs old threshold
# plot(all[1,,,1], col = c("white", "black"), breaks = c(-35, thresh$threshold, 10))

It appears that SVM finds a significant lower threshold than the amateur ‘SVM - like’ implementation (simple “best class separation”). That is while not giving the SVM the “street” - training data, because that data would be treated as an own class in SVM, which it wasn’t in the class-separation approach, where it was only regarded as “non-water” data points.

Comparison Class-Separation Threshold (left) and SVM (right)

colo = viridisLite::inferno(31)
breaks = seq(0, 31, 1)

image(class.sum, col = colo, breaks = breaks)
image(class.sum.svm, col = colo, breaks = breaks)

image(class.sum[study_area[1,]], col = colo, breaks = breaks)
image(class.sum.svm[study_area[1,]], col = colo, breaks = breaks)

plot(class.sum[moni[12,]])
plot(class.sum.svm[moni[12,]])

plot(class.sum[shape[124,]])
plot(class.sum.svm[shape[124,]])

Create Sum for Larger Area

# set Threshold -15.34187
threshold <- svmmod$x.scale$`scaled:center`
# load as shape
ext <- study_area[3,]
# bbox can be used for stars subsetting []
ext <- st_bbox(ext) # ext is xmin ymin xmax ymax
# study area one has 11860, 10.000 seems OK tile size
# calculate span
x.span <- ext[3] - ext[1] # X
y.span <- ext[4] - ext[2] # Y
# calc good number of cuts
x.cuts <- ceiling(x.span / 10000)
y.cuts <- ceiling(y.span / 10000)
# calc cut length
x.cut.by <- x.span / x.cuts
y.cut.by <- y.span / y.cuts

# tile counter
count <- 0
# go through all cuts in X direction
for (i in 1:x.cuts) {
  # go through all cuts in Y direction
  for (j in 1:y.cuts) {
    count <- count + 1
    # make extent object
    xmin <- ext[1] + (i - 1) * x.cut.by
    xmax <- ext[1] + i * x.cut.by
    ymin <- ext[2] + (j - 1) * y.cut.by
    ymax <- ext[2] + j * y.cut.by
    
    cutbox <- ext
    cutbox[1] <- xmin
    cutbox[2] <- ymin
    cutbox[3] <- xmax
    cutbox[4] <- ymax
    
    # cutbox <- extent(xmin, ymin, xmax, ymax)
    
    # load stars
    tile <- loadPolygon(cutbox)
    
    # create sum
    for (k in 1:31) {
      scene <- tile[1,,,k]
      # order is important
      scene[scene >= threshold] <- 1
      scene[scene < threshold] <- 0
      if(k == 1) {
        class.sum.svm <- scene
      }
      else {
        class.sum.svm <- class.sum.svm + scene
      }
    }
    
    # make name
    if(count < 10) {
      name <- paste0("./tiles/tile_0", count, ".tif")
    } else {
      name <- paste0("./tiles/tile_", count, ".tif")
    }
    # export
    write_stars(class.sum.svm, name)

  }
}

list <- list.files("./tiles")
for (l in 1:length(list)) {
  str <- paste0("./tiles/", list[l])
  ras <- raster(str)
  if(l < 2) {
    allRas <- ras
  }
  else {
    allRas <- raster::merge(allRas, ras)
  }
}
writeRaster(allRas, "allRas.tif", overwrite=TRUE)
allRas <- raster("allRas.tif")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved
plot(allRas)

colo = viridisLite::inferno(31)
breaks = seq(0, 31, 1)
image(allRas, col = colo, breaks = breaks)

Training of Double Bounce Class

end <- length(db$type)
for (k in 1:end) {
  if(db[k,]$type != "mixed") {
    poly <- loadPolygon(db[k,])
    # set VV, timestep
    poly <- poly[1,,,1]
    poly.df <- as.data.frame(poly)
    poly.df$type <- db[k,]$type
    poly.df <- poly.df[complete.cases(poly.df),]
    if(k == 1) {
      val.db <- poly.df[,4:5]
    }
    else {
      val.db <- rbind(val.db, poly.df[,4:5])
    }
  }
}
# boxplot
ggplot(val.db, aes(x=type, y=VV)) +
  geom_boxplot() +
  geom_hline(yintercept = -8.833493)

SVM DB Training

# prepare data
val.db$type <- as.factor(val.db$type)
attach(val.db)
## The following objects are masked from val.svm:
## 
##     type, VV
x <- subset(val.db, select=-type)
y <- type
# make model
svm.db <- svm(x,y)
summary(svm.db)
## 
## Call:
## svm.default(x = x, y = y)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1 
## 
## Number of Support Vectors:  46
## 
##  ( 28 18 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  other urban
# make prediction, confusion matrix
pred <- predict(svm.db,x)
table(pred, y)
##        y
## pred    other urban
##   other  1221     6
##   urban     2   135
# threshold is
svm.db$x.scale$`scaled:center`
## [1] -8.833493
# plot classified with threshold calculated by SVM
plot(all[1,,,1], col = c("white", "black"), breaks = c(-35, svm.db$x.scale$`scaled:center`, 20))